Imaging such as ultrasound (US) imaging has provided useful information about the interior characteristics of a subject under examination. Tissue motion estimation has been used in connection with imaging such as elastography, which has also been referred to as strain imaging, strain rate imaging, and elasticity imaging. With such imaging, stiffness or strain images of soft tissue are used to detect or classify unhealthy tissue (such as a tumor, necrotic tissue, etc.) that is stiffer relative to surrounding healthy or normal soft tissue. For this, a mechanical compression or vibration is applied to a region of interest, and the unhealthy tissue compresses less than the surrounding tissue since the strain is less than the surrounding tissue.
With axial motion (motion along the direction of the ultrasound beam), block matching, phase-shift, and frequency estimation have been used to estimate motion. Time-shift techniques, also known as speckle tracking, have used block-matching. Data used in speckle tracking can be envelope detected data (no phase information) or radio-frequency (RF) data. The RF data is data that comes before envelope-detection. For this, a one dimensional (1D) or two dimensional (2D) set of data (kernel) is chosen from a first frame. This kernel is compared to data sets of same size from a next frame within a pre-defined search region. When a best match is found, the displacement is estimated from the difference in the locations of the two regions.
Different matching algorithms have been used with block-matching. Examples of the matching algorithms have included sum of absolute differences, sum of squared differences, and maximum cross-correlation. Unfortunately, these techniques require high sampling frequency, for example, eight times the carrier frequency, which is four times the Nyquist requirement. Furthermore, the sum of absolute differences approach and the sum of squared differences approach are susceptible to identifying false minimums, and the maximum cross-correlation approach is susceptible to identifying false peaks, especially when the cross-correlation is performed on RF data.
Elastography has also required sub-sample precision of the estimates. For this, a polynomial approximation is fitted to the maximum cross-correlation or sum of absolute differences functions. One implementation of cross correlation is fastest if whole matrixes are multiplied, and the result (which is also a matrix) is compared to a temporary matrix. Adding interpolation to the cross-correlation means that the results for the different lags must be kept in memory, which increases the need for RAM. Memory access is often the bottleneck of modern systems. Using block-matching for high-precision displacement estimation requires the transfer of large amounts of data.
The phase shift estimation works exclusively with either RF data or complex-valued baseband signals, known also as in-phase and quadrature-phase signals. The fundament for phase-shift displacement estimation is the assumption that the measured signals are narrow-band and can be described as x(t)=a(t) exp(−jωt), where ω is angular frequency, to ω=2πf. The frequency of the ultrasound signal is f, and is often the carrier frequency of the pulse sent in the tissue. The time t is measured from the start of the emission for the current line. The classical estimation procedure is to take samples from the same depth from frame 1 and frame 2. It is assumed that the only difference between the two frames is due to motion.
It is also assumed that the signal x2(t) from frame 2 is a time-shifted version of the signal x1(t) from frame 1, or x2(t)=x1(t−tm). If tm is positive then the motion is away from the transducer. Assuming a narrow-band model: x1(t)=a(t)·exp(−jωt) and x2(t)=a(t−tm)·exp(−jωt+jωtm). The lag-zero cross correlation between x1(t) and x2(t) can be estimated as: R12(0)=∫a(t) exp(−jωt)·a(t)·a(tm−t) exp(jωt−jωtm) dt. The result is: R12(0)=exp(−jωm)·∫a(t)·a(tm−t). The envelope a(t) is typically a slowly varying function, and the time shift can be found from the phase of the complex correlation function R12(0): Φm=∠(R12(0))=−jωtm. The relation between depth and time is given by:
      t    =                  2        ⁢        z            c        ,where t is time, z is depth and c is speed of sound.
The displacement in depth u, which corresponds to the phase shift φm, is
  u  =            -              c                  4          ⁢          π          ⁢                                          ⁢          f                      ⁢                  ϕ        m            .      The frequency has been chosen constant, but for high-precision estimates, it must reflect the instantaneous mean frequency of the signal at a given depth. In practice, the estimation of R12(0) is done by multiplying every sample from frame 1 inside the estimation kernel with every sample from frame 2 inside the estimation kernel, located at the same place within the frame. Then, the resulting products are added together. An arctangent is used to determine the phase shift. Unfortunately, the phase shift method is susceptible to aliasing. The phase wraps around π. So a robust phase unwrapping algorithm is needed.
Another issue with this method is that the variance of the estimates increases with motion. That is, low motion is detected reliably, while large motion is not reliably detected. Variants have included zero-phase search of R12 and phase separation, where the phases of the IQ signals from the two frames are subtracted. These variants seek to overcome the aliasing problem and to decrease the variance. Their estimations start from the transducer surface, and then they keep track of the accumulated phase shift. However, if an error is made at one depth, the error may propagate to the end of the frame.
A combined autocorrelation approach estimates based on the combination of block matching and phase estimation. For both, the block matching and phase estimation, the cross-correlation function is estimated for every pixel where an estimate is needed. A correlation window from the first frame is matched against correlation windows at different displacements from the second frame. For every displacement (k, l) the normalized cross-correlation function is estimated. The maximum value of R12(k, l) is kept.
The total displacement is found from the lag at which the cross-correlation is evaluated (k, l) and from the phase of the cross-correlation at that lag:
            ϕ      =              ∠        ⁡                  (                                    R              12                        ⁡                          (                              k                ,                l                            )                                )                      ,                  ⁢                  t        m            =                        -                      ϕ                          2              ⁢              π              ⁢                                                          ⁢              f                                      +                  k                      2            ⁢                          f              s                                            ,                  ⁢    and              u      =              -                  c                      2            ⁢                          t              m                                            ,  where fs is the sampling frequency, and k is a sample in depth. This approach can be used to overcome the aliasing, as it is expected that R12 will have a maximum at a sample k within the right multiple of periods. However, the detection of the right peak of R12 requires wideband pulses, the reliable estimation of φ requires narrowband pulses, if a wrong peak is chosen, then the result shows up in the image, and this approach cannot take advantage of “envelope compression.”
Elastography is particularly sensitive to the variance of motion estimates, as strain is the first derivative of the motion. The motion within the kernel (correlation window) is not uniform. The estimation procedure aims at estimating the motion of the center of the correlation window. If there is a strong scatter at the edges of the window, then it dominates motion estimation. A periodic variation is present in the motion estimates that gives rise to an artifact known as “zebra” artifact in strain imaging. To reduce the effect of the zebra artifact, the literature suggests compressing the envelope of the RF signal prior to the estimation. The compression function known from literature is: xc(t)=sgn(x(t))×log10(1+cs×abs (x(t))), where sgn(x) returns the sign of the value and cs is a constant that controls the level of compression.